//QUIRKSIM
//Copyright (c) 2008 Gregory Andrew Robinson (garobinson@ucdavis.edu)
//
//Permission is hereby granted, free of charge, to any person
//obtaining a copy of this software and associated documentation
//files (the "Software"), to deal in the Software without
//restriction, including without limitation the rights to use,
//copy, modify, merge, publish, distribute, sublicense, and/or sell
//copies of the Software, and to permit persons to whom the
//Software is furnished to do so, subject to the following
//conditions:
//
//The above copyright notice and this permission notice shall be
//included in all copies or substantial portions of the Software.
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
//
//EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
//OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
//NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
//HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
//WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
//OTHER DEALINGS IN THE SOFTWARE.

#ifndef _PARTICLE_H
#define	_PARTICLE_H

#include <algorithm>
#include <functional>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <sys/types.h>
#include <vector>
#include "CLHEP/Vector/ThreeVector.h"

#define DISALLOW_EVIL_CONSTRUCTORS(TypeName)     \
  TypeName(const TypeName&);                     \
  TypeName& operator=(const TypeName&)
#define DISALLOW_IMPLICIT_CONSTRUCTORS(TypeName) \
  TypeName();                                    \
  DISALLOW_EVIL_CONSTRUCTORS(TypeName)

using CLHEP::Hep3Vector;

namespace quirksim {
namespace particle {

// Stores vector-valued derivatives to make integration more efficient.
struct Derivative {
  // Time derivative of momentum.
  Hep3Vector dp;
  // Time derivative of position.
  Hep3Vector dx;
};

// Base class for vector-valued forces as functions of position and velocity.
class ForceInterface {
public:
  // Parameters:
  //   x: position.
  //   v: velocity.
  virtual Hep3Vector
      GetForce(const Hep3Vector& x, const Hep3Vector& v) const = 0;
protected:
  ForceInterface() {}
  DISALLOW_EVIL_CONSTRUCTORS(ForceInterface);
};

// Base class for anything interested in watching the motion of something.
// Classes inheriting this interface should be aware of the object to which
// they are registered as a subscriber.
class ObserverInterface {
public:
  // Notify that a change occurred at a given time.
  // Parameters:
  //   time: coordinate time of the event in the observer reference frame.
  virtual void Notify(double time) = 0;
};

// Represents a single, massive, relativistic particle in 3-space.
class MassiveParticle {
public:
  // Parameters:
  //   m: velocity
  MassiveParticle(const double m);

  // Parameters:
  //   m: mass.
  //   x: position.
  //   v: velocity.
  MassiveParticle(const double m, const Hep3Vector& x, const Hep3Vector& p);

  virtual ~MassiveParticle();

  // Computes derivatives dv/dt and dx/dt.
  // Parameters:
  //   d: initial state to increment, as required by Runge-Kutta and similar.
  //   dt: time interval by which to increment.
  Derivative EvaluateDerivs(const Derivative &d, float dt) {
    Derivative output;
    output.dx = GetVelocity(momentum_ + d.dp * dt);
    output.dp = GetNetForce(position_ + d.dx * dt, output.dx);
    return output;
  }


  // Update position and momentum by numerical integration.
  // Parameters:
  //   dt: time interval.
  void Evolve(double dt);

  // Calculate total relativistic energy as a function of momentum.
  // Parameters:
  //   p: momentum used in calculation.
  double GetEnergy(const Hep3Vector& p) const;
  double GetEnergy() const {
    // E^2 = p*p*c^2 + m^2*c^4
    return GetEnergy(momentum_);
  }

  // Iterates over forces_ to find the sum.
  Hep3Vector GetNetForce(const Hep3Vector& x, const Hep3Vector& v) const;
  Hep3Vector GetNetForce() const {
    return GetNetForce(position_, GetVelocity());
  }

  // // Finds velocity as a function of momentum.
  // Parameters:
  //   p: momentum used in calculation.
  Hep3Vector GetVelocity(const Hep3Vector& p) const;
  // Finds velocity as a function of momentum.
  Hep3Vector GetVelocity() const{
    return GetVelocity(momentum_);
  }

  // Add a new force to forces_.
  void RegisterForce(ForceInterface* f);

  double mass() const { return mass_; }
  const Hep3Vector& position() const { return position_; }
  const Hep3Vector& momentum() const { return momentum_; }
  void set_position(const Hep3Vector& p) { position_ = p; }
  void set_momentum(const Hep3Vector& p) { momentum_ = p; }

protected:
  double mass_;
  Hep3Vector position_;
  Hep3Vector momentum_;
  std::vector<ForceInterface*> forces_;

private:
  DISALLOW_IMPLICIT_CONSTRUCTORS(MassiveParticle);
};

// Base class for representing the string interaction between two particles.
class QuirkString {
public:
  // Parameters:
  //   q1, q2: quirks to which the string is attached.
  QuirkString(const MassiveParticle& q1, const MassiveParticle& q2);

  // Updates position and momentum of the parameterized string.
  virtual void Evolve(double dt) = 0;

  // Returns the unit vector tangent to the string at one particle's location.
  virtual Hep3Vector GetDirection(const MassiveParticle& p) const = 0;

  // Returns the coordinate length along a string.
  virtual double GetLength() const = 0;

protected:
  const MassiveParticle& quirk1_;
  const MassiveParticle& quirk2_;
  DISALLOW_IMPLICIT_CONSTRUCTORS(QuirkString);
};

// Approximates the string interaction as a straight line.
class StraightQuirkString : public QuirkString {
public:
  StraightQuirkString(const MassiveParticle& q1, const MassiveParticle& q2);

  void Evolve(double) {}

  Hep3Vector GetDirection(const MassiveParticle& p) const {
    if (&p == &quirk1_) {
      return (quirk2_.position() - p.position()).unit();
    } else if (&p == &quirk2_) {
      return (quirk1_.position() - p.position()).unit();
    } else {
      throw std::invalid_argument("Parameter is not a particle in this string");
    }
  }

  double GetLength() const {
    return (quirk1_.position() - quirk2_.position()).mag();;
  }
};

// Bundles a pair of quirks and the string attaching them as a unit.
class QuirkPair {
public:
  // Takes particle parameters to simplify setting initial conditions.
  QuirkPair(MassiveParticle& quirk1, MassiveParticle& quirk2, double sigma);
  virtual ~QuirkPair();

  // Performs numerical integration of the system by a single step of size dt.
  void Evolve(double dt) {
    quirk1_.GetNetForce();
    quirk1_.Evolve(dt);
    quirk2_.Evolve(dt);
    string_->Evolve(dt);
  }

  // Iteratively evolves the system.
  // Parameters:
  //   dt: step size for each iteration of numerical integration.
  //   max_steps: number of integration steps to take.
  //   publish_interval: how many steps between each publication of status.
  void RunSimulation(double dt, long max_steps, long publish_interval) {
    long step;
    double time = 0;
    Publish(time);
    for(step=0; step < max_steps; ++step) {
      time += step * dt;
      Evolve(dt);
      Publish(time);
    }
  }

  Hep3Vector GetDirection(const MassiveParticle& p) const {
    Hep3Vector dir = string_->GetDirection(p);
    return dir;
  }

  // Subscribe a new observer interested in the motion of a particle.
  void Subscribe(ObserverInterface& observer);

  double sigma() const { return sigma_; }
  const MassiveParticle& quirk1() { return quirk1_; }
  const MassiveParticle& quirk2() { return quirk2_; }
  const QuirkString& string() { return *string_; }

private:
  std::vector<ObserverInterface*> observers_;
  MassiveParticle& quirk1_;
  MassiveParticle& quirk2_;
  QuirkString * const string_;
  const double sigma_;

  // Calls Notify for every observer.
  void Publish(double time) {
    std::for_each(observers_.begin(), observers_.end(),
                  std::bind2nd(std::mem_fun(&ObserverInterface::Notify), time));
  }

  DISALLOW_IMPLICIT_CONSTRUCTORS(QuirkPair);
};

// Stores the state of a QuirkPair into a QuirkPairState structure.
class QuirkPairLogger : public ObserverInterface {
public:
  // Parameters:
  //   id: arbitrary label for the output data set.
  QuirkPairLogger(std::string id, QuirkPair qp, bool print_to_screen,
                  std::ofstream pb_out);
  // Method on which we receive notification of a new piece of data.
  // Parameters:
  //   time: coordinate time of the event in the observer reference frame.
  void Notify(double time);
private:
  DISALLOW_EVIL_CONSTRUCTORS(QuirkPairLogger);
};

// Represents the Nambu-Goto string interaction for a single quirk in a pair.
class QuirkForce : public ForceInterface {
public:
  QuirkForce(const QuirkPair& qp, const MassiveParticle& p);
  Hep3Vector GetForce(const Hep3Vector& x, const Hep3Vector& v) const;
private:
  const QuirkPair& pair_;
  const MassiveParticle& quirk_;
};

} //namespace particle
} //namespace quirksim

#endif
